(******************************************************************************)
(**)                        ОТДЕЛ Собств;
(******************************************************************************
 * НАЗНАЧЕНИЕ: вычисление собственных (с.) значений и векторов
 *
 * Источники, ссылки, библиография
   Debord J., Библиотека математических подпрограмм
   tpmath.zip

   Goffe B., Программа SimAnn.FOR (Simulated Annealing in Fortran)
   http://www.netlib.org/opt/simann.f

   EISPACK: Библиотека Фортран функций для вычисления собственных значений и векторов
   http://www.netlib.org/eispack

   Marsaglia G., Тесты для генераторов случайных чисел
   http://stat.fsu.edu/~geo/diehard.html

   Moshier S., Библиотека математических подпрограмм
   http://www.moshier.net

   Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P.
   Численные рецепты на Паскале
   http://garbo.uwasa.fi/pc/turbopas.html

   Пакет численных методов для Турбо Паскаля
   Borland International, 1986
 *
 ******************************************************************************)
ИСПОЛЬЗУЕТ Матем,Вект,Матр;

ВИД 
  Вещ = Матем.Вещ;            
                                                    (* заготовки ЗАДАЧ *)
(******************************************************************************)
ЗАДАЧА^ Якоби-(A+:Матр.Вид; проходов:ЦЕЛ; точность:Вещ; СВ+:Матр.Вид; сз+:Вект.Вид):ЦЕЛ;
(* Цель:  вычисление с. значений и с. векторов методом Якоби
 ******************************************************************************
 * До:    < A >      - исходная матрица
 *        <проходов> - максимальное число проходов
 *        <точность> - требуемая точность
 * После: < СВ > - матрица нормированных с. векторов (по строкам), 
 *               у которых первый элемент > 0
 *        < сз > - с. значения
 *        < A >  - разрушенная исходная матрица
 * Ответ: 0 или НЕТ_СХОДИМОСТИ *)

(******************************************************************************)
ЗАДАЧА^ Значения-(A+:Матр.Вид; сзВещ+,сзМнм+:Вект.Вид):ЦЕЛ;
(* Цель:  вычисление с. значений квадратной матрицы
 ******************************************************************************
 * До:    < A >   - исходная матрица
 * После: <сзВещ> - вещественная часть с. значений
 *        <сзМнм> - мнимая часть с. значений
 *        < A >   - разрушенная исходная матрица
 * Ответ: 0 или НЕТ_СХОДИМОСТИ *)

(******************************************************************************)
ЗАДАЧА^ Вектор-(A-:Матр.Вид; сз,точность:Вещ; св+:Вект.Вид):ЦЕЛ;
(* Цель:  вычисление одного с. вектора для соответствующего
 *        вещественного (кратного) с. значения
 ******************************************************************************
 * До:    < A >      - исходная матрица
 *        < сз >     - с. значение
 *        <точность> - требуемая точность
 * После: < св > - нормированный с. вектор, 
 *               у которого первый элемент > 0
 * Ответ: 0 или НЕТ_СХОДИМОСТИ *)

(******************************************************************************)
ЗАДАЧА^ КорниПолинома-(Coef-:Вект.Вид; Deg:ЦЕЛ; X_Re+,X_Im+:Вект.Вид):ЦЕЛ;
(* Цель:  вычисление действительных и комплексных корней действительного
 *        полинома методом companion матрицы
 ******************************************************************************
 * До:    < Coef > - коэффициенты полинома
 *        < Deg >  - степень полинома
 * После: < X_Re > - действительные части корней (в порядке возрастания)
 *        < X_Im > - мнимые части корней
 * Ответ: 0 или НЕТ_СХОДИМОСТИ *)

(******************************************************************************)
ЗАДАЧА Якоби-(A+:Матр.Вид; проходов:ЦЕЛ; точность:Вещ; СВ+:Матр.Вид; сз+:Вект.Вид):ЦЕЛ;
(* Цель:  вычисление с. значений и с. векторов методом Якоби
 ******************************************************************************
 * До:    < A >      - исходная матрица
 *        <проходов> - наибольшее число проходов
 *        <точность> - требуемая точность
 * После: < СВ > - матрица нормированных с. векторов (по строкам), 
 *               у которых первый элемент > 0
 *        < сз > - с. значения
 *        < A >  - разрушенная исходная матрица
 * Ответ: 0 или НЕТ_СХОДИМОСТИ *)
ПЕР
  sin,cos,tg,tg2тета:Вещ;
  квCos,квSin,sincos,сумКвДиаг:Вещ;
  Aii,Ajj,Aij,Aik,Ajk,СВik,СВjk,d:Вещ;
  i,j,k,проход,посл:ЦЕЛ;
  готово:КЛЮЧ;
УКАЗ
  посл:=РАЗМЕР(A)-1;
  проход:=0;
  ОТ i:=0 ДО посл ВЫП
    ОТ j:=0 ДО посл ВЫП
      ЕСЛИ  i = j ТО
        СВ[i,j]:=1
      ИНАЧЕ
        СВ[i,j]:=0 
      КОН 
    КОН 
  КОН;

  ПОВТОРЯТЬ
    готово:=ВКЛ;
    УВЕЛИЧИТЬ(проход);
    сумКвДиаг:=0;
    ОТ i:=0 ДО посл ВЫП
      сумКвДиаг:=сумКвДиаг+Матем.кв(A[i,i])
    КОН;

    ОТ i:=0 ДО посл-1 ВЫП
      ОТ j:=i+1 ДО посл ВЫП
        ЕСЛИ  МОДУЛЬ(A[i,j]) > точность*сумКвДиаг ТО
          готово:=ОТКЛ;
          (* вычисляем поворот *)
          d:=A[i,i]-A[j,j];
          ЕСЛИ  МОДУЛЬ(d) > Матем.МАШЕПС ТО
            tg2тета:=d/(2*A[i,j]);
            tg:=Матем.знак(tg2тета)*Матем.квкор(1+Матем.кв(tg2тета))-tg2тета;
            cos:=1/Матем.квкор(1+Матем.кв(tg));
            sin:=tg*cos
          ИНАЧЕ
            cos:=Матем.КВКОР2Д2;     (* квкор(2)/2 *)
            sin:=Матем.знак(A[i,j])*Матем.КВКОР2Д2
          КОН;
          (* поворачиваем матрицу *)
          квCos:=Матем.кв(cos);
          квSin:=Матем.кв(sin);
          sincos:=sin*cos;
          Aii:=A[i,i]*квCos+2*A[i,j]*sincos+A[j,j]*квSin;
          Ajj:=A[i,i]*квSin-2*A[i,j]*sincos+A[j,j]*квCos;
          Aij:=(A[j,j]-A[i,i])*sincos+A[i,j]*(квCos-квSin);
          ОТ k:=0 ДО посл ВЫП
            ЕСЛИ  (k # i) И (k # j) ТО
              Aik:= A[i,k]*cos+A[j,k]*sin;
              Ajk:=-A[i,k]*sin+A[j,k]*cos;
              A[i,k]:=Aik;
              A[k,i]:=Aik;
              A[j,k]:=Ajk;
              A[k,j]:=Ajk
            КОН 
          КОН;
          A[i,i]:=Aii;
          A[j,j]:=Ajj;
          A[i,j]:=Aij;
          A[j,i]:=Aij;

          (* поворачиваем с. вектора *)
          ОТ k:=0 ДО посл ВЫП
            СВik:= cos*СВ[i,k]+sin*СВ[j,k];
            СВjk:=-sin*СВ[i,k]+cos*СВ[j,k];
            СВ[i,k]:=СВik;
            СВ[j,k]:=СВjk
          КОН
        КОН
      КОН
    КОН
  ДО готово ИЛИ (проход > проходов);

  (* на диагонали преобразованной матрицы находятся с. значения *)
  ОТ i:=0 ДО посл ВЫП
    сз[i]:=A[i,i]
  КОН;

  ЕСЛИ  проход > проходов ТО
    ВОЗВРАТ Матр.НЕТ_СХОДИМОСТИ
  КОН;

  (* упорядочивание с. значений и с. векторов *)
  ОТ i:=0 ДО посл-1 ВЫП
    k:=i;
    d:=сз[i];
    ОТ j:=i+1 ДО посл ВЫП
      ЕСЛИ  сз[j] > d ТО
        k:=j;
        d:=сз[j]
      КОН 
    КОН;
    Матем.обмен(сз[i],сз[k]);
    Матр.ОбменСтрок(СВ,i,k)
  КОН;

  (* первый элемент каждого с. вектора сделаем > 0 *)
  ОТ i:=0 ДО посл ВЫП
    ЕСЛИ  СВ[i,0] < 0 ТО
      ОТ j:=0 ДО посл ВЫП
        СВ[i,j]:=-СВ[i,j]
      КОН 
    КОН 
  КОН;

  ВОЗВРАТ 0
КОН Якоби;

(******************************************************************************)
ЗАДАЧА Balance(A+:Матр.Вид);
(* Цель:  сбалансировать матрицу, т.е. уменьшить ее норму, не затрагивая
 *        с. значений
 ******************************************************************************
 * До:    < A > - исходная матрица
 * После: < A > - сбалансированная матрица *)
ПОСТ
  ОСН  = 2;  (* основание счисления машинных вычислений *)
  ОСН2 = ОСН*ОСН;
ПЕР
  i,j,посл:ЦЕЛ;
  c,f,g,r,s:Вещ;
  прошли:КЛЮЧ;
УКАЗ
  посл:=РАЗМЕР(A)-1;
  ПОВТОРЯТЬ
    прошли:=ВКЛ;
    ОТ i:=0 ДО посл ВЫП
      c:=0;
      r:=0;
      ОТ j:=0 ДО посл ВЫП
        ЕСЛИ  j # i ТО
          c:=c+МОДУЛЬ(A[j,i]);
          r:=r+МОДУЛЬ(A[i,j])
        КОН 
      КОН;
      (* проверка обнуления c и r, чтобы не было переполнения *)
      ЕСЛИ  (c # 0) И (r # 0) ТО
        g:=r/ОСН;
        f:=1;
        s:=c+r;
        ПОКА c < g ВЫП
          f:=f*ОСН;
          c:=c*ОСН2
        КОН;
        g:=r*ОСН;
        ПОКА c > g ВЫП
          f:=f/ОСН;
          c:=c/ОСН2
        КОН;
        (* теперь уравниваем *)
        ЕСЛИ  (c+r)/f < 0.95*s ТО
          прошли:=ОТКЛ;
          g:=1/f;
          ОТ j:=0 ДО посл ВЫП
            A[i,j]:=A[i,j]*g 
          КОН;
          ОТ j:=0 ДО посл ВЫП
            A[j,i]:=A[j,i]*f 
          КОН
        КОН
      КОН
    КОН
  ДО прошли
КОН Balance;

(******************************************************************************)
ЗАДАЧА ElmHes(A+:Матр.Вид);
(* исключение матрицы до верхней формы Гезенберга *)
ПЕР
  i,j,m,посл:ЦЕЛ;
  x,y:Вещ;
УКАЗ
  посл:=РАЗМЕР(A)-1;
  ОТ m:=1 ДО посл-1 ВЫП
    x:=0;
    i:=m;
    ОТ j:=m ДО посл ВЫП
      ЕСЛИ  МОДУЛЬ(A[j,m-1]) > МОДУЛЬ(x) ТО
        x:=A[j,m-1];
        i:=j
      КОН 
    КОН;
    ЕСЛИ  i # m ТО
      ОТ j:=m-1 ДО посл ВЫП
        Матем.обмен(A[i,j],A[m,j]) 
      КОН;
      ОТ j:=0 ДО посл ВЫП
        Матем.обмен(A[j,i],A[j,m]) 
      КОН
    КОН;
    ЕСЛИ  x # 0 ТО
      ОТ i:=m+1 ДО посл ВЫП
        y:=A[i,m-1];
        ЕСЛИ  y # 0 ТО
          y:=y/x;
          A[i,m-1]:=y;
          ОТ j:=m ДО посл ВЫП
            A[i,j]:=A[i,j]-y*A[m,j]
          КОН;
          ОТ j:=0 ДО посл ВЫП
            A[j,m]:=A[j,m]+y*A[j,i]
          КОН
        КОН
      КОН 
    КОН
  КОН;
  ОТ i:=2 ДО посл ВЫП
    ОТ j:=0 ДО i-2 ВЫП
      A[i,j]:=0 
    КОН 
  КОН
КОН ElmHes;

(******************************************************************************)
ЗАДАЧА Hqr(A+:Матр.Вид; сзВещ+,сзМнм+:Вект.Вид):ЦЕЛ;
(* находит с. значения у верхней матрицы Гезенберга *)
ПЕР
  i,Its,j,k,l,m,n,посл:ЦЕЛ;
  Anorm,p,q,r,s,t,u,v,w,x,y,z:Вещ;

  ЗАДАЧА Sign(a,b:Вещ):Вещ;
  УКАЗ
    ЕСЛИ  b < 0 ТО ВОЗВРАТ -МОДУЛЬ(a) ИНАЧЕ ВОЗВРАТ МОДУЛЬ(a) КОН
  КОН Sign;

УКАЗ
  посл:=РАЗМЕР(A)-1;
  (* вычисление нормы *)
  k:=0;
  Anorm:=0;
  ОТ i:=0 ДО посл ВЫП
    ОТ j:=k ДО посл ВЫП
      Anorm:=Anorm+МОДУЛЬ(A[i,j]) 
    КОН;
    k:=i 
  КОН;
  n:=посл;
  t:=0;
  Its:=0;
  КОЛЬЦО
(* 60: *)
(* 70: *)
    (* поиск малого одиночного поддиагонального элемента *)
    l:=n;
    КОЛЬЦО
      ЕСЛИ  l < 1 ТО
        ВЫХОД
      КОН;
        s:=МОДУЛЬ(A[l-1,l-1])+МОДУЛЬ(A[l,l]);
      ЕСЛИ  s = 0 ТО
        s:=Anorm
      КОН;
      ЕСЛИ  s+МОДУЛЬ(A[l,l-1]) = s ТО 
        ВЫХОД
      КОН;
      УМЕНЬШИТЬ(l)
    КОН;
(* 100: *)
    ЕСЛИ  l = n ТО
(* 270: *)
      (* нашли однократный корень *)
      сзВещ[n]:=A[n,n]+t;
      сзМнм[n]:=0;
      n:=n-1;
      ЕСЛИ  n < 0 ТО
        ВОЗВРАТ 0
      КОН;
      Its:=0
      (* goto 60 *)
    ИНАЧЕ
      (* формируем сдвиг *)
      x:=A[n,n];
      y:=A[n-1,n-1];
      w:=A[n,n-1]*A[n-1,n];
      ЕСЛИ  l = n-1 ТО
(* 280: *)
        (* нашли кратный корень *)
        p:=0.5*(y-x);
        q:=Матем.кв(p)+w;
        z:=Матем.квкор(МОДУЛЬ(q));
        x:=x+t;
        ЕСЛИ  q >= 0 ТО
          (* действительноя пара *)
          z:=p+Sign(z,p);
          сзВещ[n]:=x+z;
          сзВещ[n-1]:=сзВещ[n];
          ЕСЛИ  z # 0 ТО сзВещ[n]:=x-w/z КОН;
          сзМнм[n]:=0;
          сзМнм[n-1]:=0
        ИНАЧЕ
(* 320: *)
          (* комплексная пара *)
          сзВещ[n]:=x+p;
          сзВещ[n-1]:=сзВещ[n];
          сзМнм[n]:=-z;
          сзМнм[n-1]:=z
        КОН;
(* 330: *)
        n:=n-2;
        ЕСЛИ  n < 0 ТО
          ВОЗВРАТ 0
        КОН;
        Its:=0
        (* goto 60 *)
      ИНАЧЕ
        УВЕЛИЧИТЬ(Its);
        ЕСЛИ  (Its = 10) ИЛИ (Its = 20) ТО
          (* формируем исключительный сдвиг *)
          t:=t+x;
          ОТ i:=0 ДО n ВЫП
            A[i,i]:=A[i,i]-x 
          КОН;
          s:=МОДУЛЬ(A[n,n-1])+МОДУЛЬ(A[n-1,n-2]);
          x:=0.75*s;
          y:=x;
          w:=-0.4375*Матем.кв(s)
        АЕСЛИ  Its = 30 ТО
          ВОЗВРАТ Матр.НЕТ_СХОДИМОСТИ
        КОН;
(* 130: *)
        (* поиск двух малых смежных поддиагональных элементов *)
        m:=n-2;
        КОЛЬЦО
          ЕСЛИ  m < l ТО
            ВЫХОД
          КОН;
          z:=A[m,m];
          r:=x-z;
          s:=y-z;
          p:=(r*s-w)/A[m+1,m]+A[m,m+1];
          q:=A[m+1,m+1]-z-r-s;
          r:=A[m+2,m+1];
          s:=МОДУЛЬ(p)+МОДУЛЬ(q)+МОДУЛЬ(r);
          p:=p/s;
          q:=q/s;
          r:=r/s;
          ЕСЛИ  m = l ТО 
            ВЫХОД
          КОН;
          v:=МОДУЛЬ(p)*(МОДУЛЬ(A[m-1,m-1])+МОДУЛЬ(z)+МОДУЛЬ(A[m+1,m+1]));
          u:=МОДУЛЬ(A[m,m-1])*(МОДУЛЬ(q)+МОДУЛЬ(r));
          ЕСЛИ  u+v=v ТО 
            ВЫХОД
          КОН;
          УМЕНЬШИТЬ(m)
        КОН;
(* 150: *)
        ОТ i:=m+2 ДО n ВЫП
          A[i,i-2]:=0;
          ЕСЛИ  i # (m+2) ТО 
            A[i,i-3]:=0
          КОН 
        КОН;
        (* двойной QR шаг разрешая строки l на n и столбцы m на n *) 
        ОТ k:=m ДО n-1 ВЫП
          ЕСЛИ  k # m ТО
            p:=A[k,k-1];
            q:=A[k+1,k-1];
            r:=0;
            ЕСЛИ  k # (n-1) ТО
              r:=A[k+2,k-1]
            КОН;
            x:=МОДУЛЬ(p)+МОДУЛЬ(q)+МОДУЛЬ(r);
            ЕСЛИ  x # 0 ТО
              p:=p/x;
              q:=q/x;
              r:=r/x
            КОН 
          КОН;
(* 170: *)
          s:=Sign(Матем.квкор(Матем.кв(p)+Матем.кв(q)+Матем.кв(r)),p);
          ЕСЛИ  s # 0 ТО
            ЕСЛИ  k = m ТО
              ЕСЛИ  l # m ТО
                A[k,k-1]:=-A[k,k-1]
              КОН
            ИНАЧЕ
              A[k,k-1]:=-s*x
            КОН;
            p:=p+s;
            x:=p/s;
            y:=q/s;
            z:=r/s;
            q:=q/p;
            r:=r/p;
            (* изменение строк *)
            ОТ j:=k ДО n ВЫП
              p:=A[k,j]+q*A[k+1,j];
              ЕСЛИ  k # (n-1) ТО
                p:=p+r*A[k+2,j];
                A[k+2,j]:=A[k+2,j]-p*z
              КОН;
              A[k+1,j]:=A[k+1,j]-p*y;
              A[k,j]:=A[k,j]-p*x
            КОН;
            (* изменение столбцов *)
            ОТ i:=l ДО Матем.минЦ(n,k+3) ВЫП
              p:=x*A[i,k]+y*A[i,k+1];
              ЕСЛИ  k # (n-1) ТО
                p:=p+z*A[i,k+2];
                A[i,k+2]:=A[i,k+2]-p*r
              КОН;
              A[i,k+1]:=A[i,k+1]-p*q;
              A[i,k]:=A[i,k]-p
            КОН 
          КОН 
        КОН
      КОН 
    КОН 
  КОН
КОН Hqr;

(******************************************************************************)
ЗАДАЧА Значения-(A+:Матр.Вид; сзВещ+,сзМнм+:Вект.Вид):ЦЕЛ;
(* Цель:  вычисление с. значений квадратной матрицы
 ******************************************************************************
 * До:    < A >   - исходная матрица
 * После: <сзВещ> - вещественная часть с. значений
 *        <сзМнм> - мнимая часть с. значений
 *        < A >   - разрушенная исходная матрица
 * Ответ: 0 или НЕТ_СХОДИМОСТИ *)
УКАЗ
  Balance(A);
  ElmHes(A);
  ВОЗВРАТ Hqr(A,сзВещ,сзМнм)
КОН Значения;

(******************************************************************************)
ЗАДАЧА Вектор-(A-:Матр.Вид; сз,точность:Вещ; св+:Вект.Вид):ЦЕЛ;
(* Цель:  вычисление одного с. вектора для соответствующего
 *        вещественного (кратного) с. значения
 ******************************************************************************
 * До:    < A >    - исходная матрица
 *        < сз >   - с. значение
 *        <точность> - требуемая точность
 * После: < св >     - нормированный с. вектор, 
 *                   у которого первый элемент > 0
 * Ответ: 0 или НЕТ_СХОДИМОСТИ *)
ПЕР
  код,i,посл:ЦЕЛ;
  A1:Матр.Доступ;

  ЗАДАЧА SetMatrix(A-,A1+:Матр.Вид; сз:Вещ);
  (* строит A1 = A-сз*i *)
  ПЕР
    i:ЦЕЛ;
  УКАЗ
    Матр.Списать(A,A1);
    ОТ i:=0 ДО РАЗМЕР(A)-1 ВЫП
      A1[i,i]:=A[i,i]-сз 
    КОН
  КОН SetMatrix;

  ЗАДАЧА Solve(A-:Матр.Вид; n:ЦЕЛ; точность:Вещ; св+:Вект.Вид):ЦЕЛ;
  (* решает систему A*x = 0 после закрепления за n-м неизвестным 1 *)
  ПЕР
    A1,W:Матр.Доступ;
    b,s,x:Вект.Доступ;
    код,i,i1,j,j1,посл:ЦЕЛ;
  УКАЗ
    посл:=РАЗМЕР(A)-1;
    СОЗДАТЬ(A1,посл,посл);
    СОЗДАТЬ(W,посл,посл);
    СОЗДАТЬ(b,посл);
    СОЗДАТЬ(s,посл);
    СОЗДАТЬ(x,посл);

    i1:=0;
    ОТ i:=0 ДО посл ВЫП
      ЕСЛИ  i # n ТО
        j1:=0;
        ОТ j:=0 ДО посл ВЫП
          ЕСЛИ  j # n ТО
            A1[i1,j1]:=A[i,j];
            УВЕЛИЧИТЬ(j1)
          КОН
        КОН;
        b[i1]:=-A[i,n];
        УВЕЛИЧИТЬ(i1)
      КОН 
    КОН;

    код:=Матр.РазложитьНаSV(A1^,s^,W^);

    ЕСЛИ  код = 0 ТО
      Матр.ОбнулитьS(s^,точность);
      Матр.РешитьИзSV(A1^,s^,W^,b^,x^);

      (* обновить с. вектор *)
      i1:=0;
      ОТ i:=0 ДО посл ВЫП
        ЕСЛИ  i = n ТО
          св[i]:=1
        ИНАЧЕ
          св[i]:=x[i1];
          УВЕЛИЧИТЬ(i1)
        КОН 
      КОН
    КОН;
    ВОЗВРАТ код;
  КОН Solve;

  ЗАДАЧА ZeroVector(b-:Вект.Вид; точность:Вещ):КЛЮЧ;
  (* ВКЛ, если вектор b 0-й *)
  ПЕР
    i:ЦЕЛ;
  УКАЗ
    ОТ i:=0 ДО РАЗМЕР(b)-1 ВЫП
      ЕСЛИ  МОДУЛЬ(b[i]) >= точность ТО
        ВОЗВРАТ ОТКЛ
      КОН
    КОН;
    ВОЗВРАТ ВКЛ
  КОН ZeroVector;

  ЗАДАЧА CheckEigenVector(A1-:Матр.Вид; св-:Вект.Вид; точность:Вещ):КЛЮЧ;
  (* ВКЛ, если выполняется равенство A1*св = 0 *)
  ПЕР
    i,k,посл:ЦЕЛ;
    b:Вект.Доступ;
  УКАЗ
    посл:=РАЗМЕР(A1)-1;
    СОЗДАТЬ(b,посл+1);
    (* строим b = A1*св *)
    ОТ i:=0 ДО посл ВЫП
      ОТ k:=0 ДО посл ВЫП
        b[i]:=b[i]+A1[i,k]*св[k]
      КОН 
    КОН;
    ВОЗВРАТ ZeroVector(b^,точность)
  КОН CheckEigenVector;

  ЗАДАЧА Normalize(св+:Вект.Вид);
  (* нормируем с. вектор и устанавливаем первый элемент >= 0 *)
  ПЕР
    сумма,Norm:Вещ;
    i,посл:ЦЕЛ;
  УКАЗ
    посл:=РАЗМЕР(св)-1;
    сумма:=0;
    ОТ i:=0 ДО посл ВЫП
      сумма:=сумма+Матем.кв(св[i]) 
    КОН;
    Norm:=Матем.квкор(сумма);
    ОТ i:=0 ДО посл ВЫП
      св[i]:=св[i]/Norm 
    КОН;
    ЕСЛИ  св[0] < 0 ТО
      ОТ i:=0 ДО посл ВЫП
        св[i]:=-св[i]
      КОН 
    КОН
  КОН Normalize;

УКАЗ
  посл:=РАЗМЕР(A)-1;
  СОЗДАТЬ(A1,посл+1,посл+1);

  (* строим A1 = A-сз*i *)
  SetMatrix(A,A1^,сз);

  (* попытаться решить систему A1*св=0 исключая 1 уравнение *)
  ОТ i:=0 ДО посл ВЫП
    ЕСЛИ  (Solve(A1^,i,точность,св) = 0) И CheckEigenVector(A1^,св,точность) ТО
      Normalize(св);
      ВОЗВРАТ 0
    КОН
  КОН;
  ВОЗВРАТ Матр.НЕТ_СХОДИМОСТИ 
КОН Вектор;

(******************************************************************************)
ЗАДАЧА КорниПолинома-(Coef-:Вект.Вид; Deg:ЦЕЛ; X_Re+,X_Im+:Вект.Вид):ЦЕЛ;
(* Цель:  вычисление действительных и комплексных корней действительного
 *        полинома методом companion матрицы
 ******************************************************************************
 * До:    < Coef > - коэффициенты полинома
 *        < Deg >  - степень полинома
 * После: < X_Re > - действительные части корней (в порядке возрастания)
 *        < X_Im > - мнимые части корней
 * Ответ: 0 или НЕТ_СХОДИМОСТИ *)
ПЕР
  A:Матр.Доступ;   (* companion матрица *)
  n:ЦЕЛ;           (* размер матрицы  *)
  i,j,k:ЦЕЛ;  
  код:ЦЕЛ;
  врем:Вещ;
УКАЗ
  n:=Deg-1;
  СОЗДАТЬ(A,n+1,n+1);

  (* заполняем companion матрицу *)
  ОТ j:=0 ДО n ВЫП
    A[0,j]:=-Coef[Deg-j-2]/Coef[Deg-1]
  КОН;
  ОТ j:=0 ДО n-1 ВЫП
    A[j+1,j]:=1 
  КОН;

  (* корни полинома являются и с. значениями companion матрицы *)
  Balance(A^);
  код:=Hqr(A^,X_Re,X_Im);

  ЕСЛИ  код = 0 ТО
    (* упорядочивание корней в порядке возрастания действительных частей *)
    ОТ i:=0 ДО n-1 ВЫП
      k:=i;
      врем:=X_Re[i];
      ОТ j:=i+1 ДО n ВЫП
        ЕСЛИ  X_Re[j] < врем ТО
          k:=j;
          врем:=X_Re[j]
        КОН 
      КОН;
      Матем.обмен(X_Re[i],X_Re[k]);
      Матем.обмен(X_Im[i],X_Im[k])
    КОН
  КОН;
  ВОЗВРАТ код
КОН КорниПолинома;

КОН Собств.